Self-structuring of Granular Media under Internal Avalanching 
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We study the phenomenon of internal avalanching within the context of recently proposed "Tetris" 
lattice models for granular media. We define a recycling dynamics under which the system reaches 
a steady state which is self-structured, i.e. it shows a complex interplay between textured internal 
structures and critical avalanche behavior. Furthermore we develop a general mean-field theory for 
this class of systems and discuss possible scenarios for the breakdown of universality. 
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There has been a lot of interest in understanding the 
internal structure and geometry of granular packings 
0. The rich phenomenology observed in experiments 
studying compaction, segregation and force distributions, 
amongst other things, has prompted a number of numeri- 
cal and analytical studies. In another context the interest 
in granular media has also been triggered by the search 
for self-organised criticality (SOC) ||. As a result sur- 
face avalanches in granular piles have been extensively 
studied both experimentally Q and in computer models 
in order to identify clear regimes of SOC-like be- 
havior ||. In this Letter, we focus our interest on the 
interplay between the internal structure of packings and 
power-law avalanche distributions. We define a steady 
state dynamics under which the medium reaches a self- 
structured critical state. We also focus on the internal 
structure of this state and find, very interestingly, that 
it can be highly inhomogeneous with strong segregation 
and ordering effects. 

The model we have investigated is a recently proposed 
simple lattice model for describing slow dynamical pro- 
cesses in granular media . The basic ingredients of this 
model are the geometric constraints involved in packing 
particles of different shapes. This model is seen to repro- 
duce experimentally observed phenomena such as slow 
relaxation in compaction segregation Q|, as well as 
aging @. 

Within the context of this model, we study the phe- 
nomenon of internal avalanches occurring under small 
perturbations. But as opposed to previous works 
we focus on the stationary state that a system reaches 
under the continued process of removing a particle from 
the bottom layer and adding it back to the top of the 
system. Under this dynamics the system reaches a well 
defined "critical" steady state in which the avalanche dis- 
tributions decay as power laws. Most interestingly, we 
find that in order to achieve this effect, the system re- 
structures under this dynamics to a very inhomogeneous 



state with ordered regions (grains) separated by disor- 
dered low-density channels (grain boundaries) which act 
as preferential pathways for these avalanches. We per- 
form our numerical experiments for particles of several 
different shapes and find that the steady state reached 
is always as described above, with an exponent for the 
power-law distribution which is the same for a large class 
of particle shapes. We furthermore develop a mean-field 
theory for systems undergoing this dynamics and explain, 
within this context, why we observe a universal power- 
law distribution. We elaborate on this point by consider- 
ing a case when an important change in the rules of sta- 
bility of particles changes the steady state reached and 
hence the universality class of the phenomenon. 

We briefly review the definitions and some basic prop- 
erties of the Tetris models || used in our simulations. 
Frustration arises in granular packings owing to excluded 
volume effects of particles of different shapes. This ge- 
ometrical feature is captured in the Tetris model. In 
the following, we present results for the simplest version 
of the model where the particles are either rods with 
two kinds of orientations, more complicated shapes such 
as "T"- shaped particles with two kinds of orientations 
or "crosses" with arms of randomly distributed lengths 
in the framework of the so-called Random Tetris Model 
(RTM) go). 

The Tetris model can be defined as a system of parti- 
cles which occupy the sites of a square lattice tilted by 
45° with periodic boundary conditions in the horizontal 
direction (cylindrical geometry) and a rigid wall at the 
bottom. Particles cannot overlap and this condition pro- 
duces very strong constraints (frustration) on their rela- 
tive positions. This is illustrated for "T-shaped" particles 
in Fig. ([!]). In general each particle can be schematized as 
a cross with arms of different lengths which can be cho- 
sen in a regular fla] or in a random way JhJ. The system 
is initialized by inserting the particles at the top of the 
system, one at a time, and letting them move down under 
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gravity. The particles perform an oriented random walk 
on the lattice until they reach a stable position defined 
as a position from which they cannot fall any further be- 
cause of other particles below them. The particles retain 
their orientations as they move, i.e. they are not allowed 
to rotate. We now introduce the following dynamics un- 
der which the system evolves. A particle is removed from 
a random position at the base. This could destabilise its 
neighbouring particles above one of which may then fall 
down if the geometry of the packing allows for the mo- 
tion (i.e., if the orientation of the particle fits the local 
conformation). In this case, the disturbance propagates 
upwards destabilizing particles in the layer above and 
so on. We update the system sequentially, moving all 
the unstable particles down till the system is once more 
stable. The removed particle is then added back at a 
random position at the top of the system. This process 
is continued till the system reaches a steady state. 

Similar procedures have been studied before for other 
models pJE^. While long-ranged avalanche distribu- 
tions have been found in |llf| , the update rule assumed in 
Jl2| does not lead to a critical state. We go beyond these 
previous works by studying here, in detail, the interplay 
between the avalanche distribution and the density pro- 
file of the medium. We explain the means by which the 
system reaches a critical state by developing a generic 
mean-field theory for avalanche distributions in dense or 
loose packings. We utilise the possibility afforded by this 
model, of easily changing particle shapes, to study this 
behaviour for a wide variety of particle shapes. Most in- 
terestingly we also find that this "critical" steady state 
is inhomogeneous and strongly ordered, different from 
those ordinarily studied in most SOC systems. These 
are thus some of the new features reported in the present 
study. 

Fig. 2a is a picture of a packing of in the steady state. 
As can be seen it shows a complex textured structure. 
Namely, beginning from an initial state in which particles 
of different shapes are homogeneously mixed, the packing 
always "segregates" under the dynamics so as to form or- 
dered high density grains separated by grain boundaries 
at lower densities. All avalanches preferentially propa- 
gate inside these grain boundaries i.e. no matter where 
the initial seed, the avalanches find their way into the 
boundary region (see Fig-d-left)). 

The size of an avalanche is defined as the total num- 
ber of particles destabilized by the process of removing 
a particle at the bottom. The size distribution of the 
avalanches decays like a power P(s) ~ s~ T . This was 
studied for the three different types of particles described 
above. Time averages were performed in the steady state 
over ~ 10 6 configurations in order to obtain good statis- 
tics. Fig.(||) shows the avalanche distributions obtained 
for two different choices of the particles: the T's shown 
in Fig.(^) and particles with random shapes obtained in 
the framework of the RTM. In both cases one observes a 



scaling behavior for the avalanches with t — 1.5 ± 0.05. 
In the case of the rods, the result is sensitive to the aspect 
ratio of the system but for systems of about equal width 
and height the exponent of the avalanche distribution is 
again the same. 

We now turn to a discussion of the avalanche statistics 
within the framework of a mean-field theory that we de- 
velop for this class of systems. It is apparent from our 
numerical studies reported so far that under this dynam- 
ics, the system reaches a steady state which is critical. 
The reason could be the following. Taking out a particle 
in the last layer creates a void in the packing. This can ei- 
ther move up (by exchanging place with a particle), die (if 
nothing above is destabilized) or free another neighboring 
void (and hence multiply) and propagate. The dynamics 
is thus essentially like a branching-annihilating process 
on the lattice where the probabilities for annihilating Pq, 
branching P2 and propagation Pi = I — P0—P2 depend on 
the density of the packing. However there is also a feed- 
back effect. The avalanche distribution can in its turn 
affect the density of the system : large avalanches that 
reach the top tend to compactify the system and small 
avalanches make the system looser. 

We can make the above arguments more precise in the 
following way. Let ph(t) be the cumulative density of the 
system up to height h at time t. Then the density of the 
system at time t + 1 will be: 

Ph(t + 1) - Ph (t) = -l/L 2 + (t) * h 1 / 1? (1) 

where L is the linear size (height) of the system. The 
first term of the RHS represents the effect of removing 
one particle. This is the sole contribution of avalanches 
which die before reaching the height h. On the other 
hand, those avalanches which reach at least a height h 
have the additional effect of increasing the density of the 
system by an amount equal to the number of voids which 
escape at h. This is equal to the width of the avalanche at 
h. If the avalanches are self-affine (as in pij] , and also in 
the case studied here), i.e. an avalanche of height h has 
a width of h 1 , then the density increase is precisely given 
by the second term. The coefficient a(t) is just a random 
variable which takes the value 1 whenever an avalanche 
reaches at least a height h and a value otherwise. In the 
steady state, we can perform a time average on Eq. ([[]). 
We expect the LHS to vanish in this case. To evaluate 
the RHS, we note that the time average of a(t) is simply 
1 — P(h)dh. We measure P(h) numerically and ob- 
serve the existence of a scaling region where P(h) ~ 
with f3 — 1.95 ± 0.05. We also independently measure 
7 (by measuring the characteristic size s* — h 1+1 of the 
avalanche size cut-off at height h). We find 7 = 0.9 ±0.1. 

The above equation makes a prediction for the 
avalanche exponent. The steady state condition requir- 
ing that the average density of the system (p) = const 
implies that — 1 =7- From the numerically mea- 
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sured value of mentioned above, we see that we find 
7 ~ 0.95 consistent with the numerically measured value 
of 7. Making a change of variables from the avalanche 
height h to the avalanche size s gives us the relation 
(1 -P {h))li< = /iV^ (1+7)(1 ~ r) = 1 where P(s) ~ s~ T is 
the avalanche size distribution in the steady state. The 
above scaling relation for j3 — 1 = 7 then translates to 
(also obtained in |ll|] using a steady state argument) 

t = 1 + 7/(1+7). (2) 

Using again our numerical estimate for 7 we find r ~ 
1.47 ± 0.05 consistent with the data shown in Fig. ||. 

A more complete and self-consistent description of the 
observed phenomenology can be obtained complementing 
equation (|l|) with an equation for the avalanche distribu- 
tion P(s) in terms of p, i.e. with an equation 

P(s(t)) = F(p(t)) (3) 

where F indicates a generic function of p{t). The two 
coupled equations ([!]) and (||) should then describe the 
evolution of the system to a steady state given by a criti- 
cal density p c with an avalanche distribution decaying as 
a power-law. In general, it is difficult to write an exact 
equation for the avalanche distribution in terms of the 
density except for avalanches propagating on the Bethe 
lattice In this case, it is possible to show quite sim- 
ply that the feedback effect of Eq. (||) on Eq. (|l|) results 
in the system reaching a critical density p c given sim- 
ply by the equation P\(p c ) + 2 * P2(p c ) = 1 where Pi 
and Pi are the probabilities for propagation and branch- 
ing respectively, introduced before. We have investigated 
analytically and numerically that the mean-field theory 
is insensitive to the exact functional form of the birth- 
death probabilities and avalanches always decay with an 
exponent r = 1.5 at the critical density. A more detailed 
analysis of the above equations considering different ex- 
plicit forms of F is considered elsewhere fl4fl . 

The scaling relation (||) always holds for systems with 
open boundary conditions provided that there is a com- 
pact bulk packing. This poses an upper limit to the ex- 
ponent t. For non-fractal bulk packings with a smooth 
free surface, 7 cannot be larger than 1 and hence r can- 
not be larger than 1.5. It is interesting to note that the 
avalanches decay with the same exponent as in mean- 
field theory. However, the reason for this exponent here 
is that the avalanches propagate in a conical region (im- 
plied by 7 ~ 1 ) centered around the grain boundary 
(since, as mentioned before, avalanches propagate most 
of the time at the grain boundary). These facts imply, 
from the scaling relation (|J), that r = 1.5. 

Although our results have so far shown a universal be- 
havior, we identify within the framework of this theory 
at least one clear instance of the breakdown of this uni- 
versality. This has to do with having a very loose pack- 
ing in the system. If this is the case, particles can fall 



large distances in the course of an avalanche and com- 
pactify the system far below. It would then not only be 
the width of the avalanche at height h which would con- 
tribute to the compactification but some fraction of the 
whole avalanche above h. Such an effect is clearly not 
taken into account in Eq.(Q) which hence implicitly as- 
sumes that particles only fall short distances. We thus 
have to rewrite the above mean-field theory for a loose 
system for which the particles can fall large distances. We 
can quantify the above statements by rewriting Eq. (|l|) in 
the following manner: 

Ph(t + 1) - Ph {t) = + Jjs ~ s*rs- T ds (4) 

where s* is the typical size of an avalanche reaching a 
height h and a is a measure of how much of this avalanche 
contributes to heights less than h. Making a change of 
variables and taking the s* dependence out of the in- 
tegral, we find that the relevant scaling relation is now 
expressed in terms of s* as t = 1 + a. Since a < 1 (if 
the total avalanche above h contributes then a = 1) we 
find that for systems with very loose packings the upper 
bound for the avalanche distribution is now r < 2 and 
not 1.5 as before. 

We have checked this by changing the stability condi- 
tion for particles to get a much looser packing. In all the 
cases considered above, the particles need to be stable 
in two directions in order not to fall. We modified this 
by looking at a system of sticky particles, in which one 
downward contact in either direction suffices for stability. 
Repeating the same recycling procedure used throughout 
the paper we find in this case a stationary state with a 
non-compact bulk packing (Fig.^-right). For this system, 
we find an avalanche distribution in the steady state char- 
acterized by an exponent r = 1.9 + 0.05 (see Fig]|) out of 
the range of validity of Eq. (§) and in the range of validity 
of the scaling relation predicted by Eq. (H|). 

There are several features that it is of interest to inves- 
tigate further. An instability mechanism for producing 
structured steady states has yet to be developed [fl4"|| . 
Further it would be interesting to see how these struc- 
tures coexist with power-law avalanches and whether fi- 
nite driving destroys this effect. In this context it can be 
seen from Fig. || that the big avalanches are enhanced 
well over the power-law. It could be of interest to inves- 
tigate whether this is just a finite size effect or whether 
the structures play a role in this |l^] . Within the context 
of this model we have also studied a system of spheri- 
cal particles (i.e. crosses with roughly equal extensions 
in either direction). This is the case closest to the one 
studied in |ll[ in which the particles are all the same 
shape. We find that though a density plot is not suffi- 
cient to spot structures, an activity plot (marking how 
many avalanches pass through every site over a period of 
time) shows very distinctly that there are always long- 
lived loose regions where avalanches preferentially prop- 
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agate with r ~ 1.5. It is hence tempting to conclude 
that this dynamics always results in long-lived inhomo- 
geneities (with easy channels for particle flow) which af- 
fect the avalanche distribution. Finally, it is interesting 
to speculate what our results imply for possible exper- 
iments on the phenomenon of internal avalanches. One 
implication might be that a real system subjected to the 
continual process of removal and addition of grains will 
"fracture" (as in our model) developing easy regions for 
particle flow. It would be very interesting to see whether 
this is observable. 

Acknowledgements: HJH, SK, SSM and SR would 
like to thank CEFIPRA. In particular, SK and SSM 
would like to acknowledge financial support from CE- 
FIPRA under project no 1508-3/192. VL acknowledges 
financial support under project ERBFMBICT961220. 
This work has also been partially supported from the 
European Network-Fractals under contract No. FM- 
RXCT980183. 



[1] For a recent introduction to the overall phenomenology 
see Proceedings of the NATO Advanced Study Institute 
on Physics of Dry Granular Media, Eds. H. J. Herrmann 
et al, Kluwer Academic Publishers, Netherlands (1998). 

[2] P. Bak, C. Tang and K. Weisenfeld, Phys. Rev. Lett. 59, 
381 (1987). 

[3] H. M. Jaeger, C. -H. Liu and S. R. Nagel, Phys. Rev. Lett. 

62, 40 (1989); G. A. Held et al, Phys. Rev. Lett. 65, 1120 

(1990); M. Bretz et al, Phys. Rev. Lett. 69, 2431 (1992). 
[4] D. Dhar, Phys. Rev. Lett. 64, 1613 (1990); L. P. Kadanoff 

et al, Phys. Rev . A. 39, 6524 (1989); S. S. Manna, J. 

Phys. A 24, L363 (1992). 
[5] V. Frette, K. Christensen, A. Malthe-S0renssen, J. Feder, 

T. J0ssang and P. Meakin, Nature 379, 49 (1996). 
[6] E. Caglioti, V. Loreto, H.J. Herrmann and M. Nicodemi, 

Phys. Rev. Lett. 79, 1575 (1997). 
[7] E. Caglioti, A. Coniglio, H.J. Herrmann, V. Loreto and 

M. Nicodemi, Europhys. Lett. 43, 591 (1998). 
[8] M. Nicodemi and A. Coniglio, Phys. Rev. Lett. 82, 916 

(1999). 

[9] S. Krishnamurthy, H. J. Herrmann, V. Loreto, M. 

Nicodemi and S. Roux, to appear in Fractals (1998). 
[10] E. Caglioti, S. Krishnamurthy and V. Loreto, Random 

Tetris Model for Granular Systems, unpublished (1998). 
[11] R. E. Snyder and R. C. Ball, Phys. Rev. E, 49, 104 

(1994). 

[12] S. S. Manna, Proceedings of the International Confer- 
ence on Statistical Physics, Calcutta, 1999, to appear in 
Physica A. 

[13] This mean-field is similar in spirit to the self-organised 
branching process studied in Zapperi et al., Phys. Rev. 
Lett. 75, 4071 (1995). 

[14] S. Krishnamurthy, V. Loreto and S. Roux, unpublished 
(1999). 




FIG. 1. A steady state configuration of the Tetris model 
with "T-shaped" particles with two different orientations. 
The boundary conditions are periodic in the horizontal di- 
rection. The black particles are those which rearrange in the 
avalanche caused by removing the lowest particle. 
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FIG. 2. Typical avalanches in the steady state for a system 
of T-shaped particles (left) and sticky particles (right). 
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FIG. 3. P(s) vs s in the steady state for "T-shaped" par- 
ticles and the "crosses" (RTM): r = 1.5 ± 0.05 in both 
cases. The system sizes shown are Lx = 100, Ly — 500 and 
Lx = 200, Ly = 650, 1000 respectively for the "Tees" and 
Lx = 100, Ly = 150 , Lx = 200, Ly = 300 for the "crosses". 
The last curve shows the avalanche distribution for a system of 
sticky particles (see text). In this case one gets r = 1.9±0.05. 
The system size shown is Lx = 200, Ly = 350. 
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